# Third tutorial¶

## Crystalline silicon.¶

This tutorial aims at showing you how to get the following physical properties, for an insulator:

• the total energy
• the lattice parameter
• the band structure (actually, the Kohn-Sham band structure)

You will learn about the use of k-points, as well as the smearing of the plane-wave kinetic energy cut-off.

Visualisation tools are NOT covered in this tutorial. Powerful visualisation procedures have been developed in the Abipy context, relying on matplotlib. See the README of Abipy and the Abipy tutorials.

This tutorial should take about 1 hour.

Note

Supposing you made your own install of ABINIT, the input files to run the examples are in the ~abinit/tests/ directory where ~abinit is the absolute path of the abinit top-level directory. If you have NOT made your own install, ask your system administrator where to find the package, especially the executable and test files.

To execute the tutorials, create a working directory (Work*) and copy there the input files and the files file of the lesson. This will be explicitly mentioned in the first lessons, that will tell you more about the files file (see also section 1.1). The files file ending with _x (e.g. tbase1_x.files) must be edited every time you start to use a new input file.

Most of the tutorials do not rely on parallelism (except specific tutorials on parallelism). However you can run most of the tutorial examples in parallel, see the topic on parallelism.

In case you work on your own PC or workstation, to make things easier, we suggest you define some handy environment variables by executing the following lines in the terminal:

export ABI_HOME=Replace_with_the_absolute_path_to_the_abinit_top_level_dir
export PATH=$ABI_HOME/src/98_main/:$PATH
export ABI_TESTS=$ABI_HOME/tests/ export ABI_PSPDIR=$ABI_TESTS/Psps_for_tests/  # Pseudopotentials used in examples.


Examples in this tutorial use these shell variables: copy and paste the code snippets into the terminal (remember to set ABI_HOME first!). The ‘export PATH’ line adds the directory containing the executables to your PATH so that you can invoke the code by simply typing abinit in the terminal instead of providing the absolute path.

## Computing the total energy of silicon at a fixed number of k-points¶

Before beginning, you might consider working in a different subdirectory, as for tutorial 1 or 2. Why not Work3?

The file tbase3_x.files lists the file names and root names. You can copy it in the Work3 directory and change it as you did in the first and second tutorials. You can also copy the file tbase3_1.in inside the Work3 directory with:

while $ABI_TESTS/tutorial/Refs/tbase3_5.out is a reference output file. You should find the band structure starting at (second dataset):  Eigenvalues ( eV ) for nkpt= 40 k points: kpt# 1, nband= 8, wtk= 1.00000, kpt= 0.5000 0.0000 0.0000 (reduced coord) -3.78815 -1.15872 4.69668 4.69668 7.38795 9.23867 9.23867 13.45707 kpt# 2, nband= 8, wtk= 1.00000, kpt= 0.4500 0.0000 0.0000 (reduced coord) -3.92759 -0.95774 4.71292 4.71292 7.40692 9.25561 9.25561 13.48927 kpt# 3, nband= 8, wtk= 1.00000, kpt= 0.4000 0.0000 0.0000 (reduced coord) -4.25432 -0.44393 4.76726 4.76726 7.46846 9.31193 9.31193 13.57737 kpt# 4, nband= 8, wtk= 1.00000, kpt= 0.3500 0.0000 0.0000 (reduced coord) -4.64019 0.24941 4.85732 4.85732 7.56855 9.38323 9.38323 13.64601 ....  One needs a graphical tool to represent all these data. In a separate file (_EIG), you will find the list of k-points and the eigenenergies (the input variable prteig is set by default to 1). Even without a graphical tool we will have a quick look at the values at L, $\Gamma$, X and $\Gamma$ again:  kpt# 1, nband= 8, wtk= 1.00000, kpt= 0.5000 0.0000 0.0000 (reduced coord) -3.78815 -1.15872 4.69668 4.69668 7.38795 9.23867 9.23867 13.45707 kpt# 11, nband= 8, wtk= 1.00000, kpt= 0.0000 0.0000 0.0000 (reduced coord) -6.17005 5.91814 5.91814 5.91814 8.44836 8.44836 8.44836 9.17755 kpt# 23, nband= 8, wtk= 1.00000, kpt= 0.0000 0.5000 0.5000 (reduced coord) -1.96393 -1.96393 3.00569 3.00569 6.51173 6.51173 15.95524 15.95524 kpt# 39, nband= 8, wtk= 1.00000, kpt= 1.0000 1.0000 1.0000 (reduced coord) -6.17005 5.91814 5.91814 5.91814 8.44836 8.44836 8.44836 9.17755  The last $\Gamma$ is exactly equivalent to the first $\Gamma$. It can be checked that the top of the valence band is obtained at $\Gamma$ (=5.91814 eV). The width of the valence band is 12.09 eV, the lowest unoccupied state at X is 0.594 eV higher than the top of the valence band, at $\Gamma$. The Si is described as an indirect band gap material (this is correct), with a band-gap of about 0.594 eV (this is quantitatively quite wrong: the experimental value 1.17 eV is at 25 degree Celsius). The minimum of the conduction band is even slightly displaced with respect to X, see kpt # 21. This underestimation of the band gap is well-known (the famous DFT band-gap problem). In order to obtain correct band gaps, you need to go beyond the Kohn-Sham Density Functional Theory: use the GW approximation. This is described in the first tutorial on GW. For experimental data and band structure representation, see the book by M.L. Cohen and J.R. Chelikowski [Cohen1988]. Important There is a subtlety that is worth to comment about. In non-self-consistent calculations, like those performed in the present band structure calculation, with iscf = -2, not all bands are converged within the tolerance tolwfr. Indeed, the two upper bands (by default) have not been taken into account to apply this convergence criterion: they constitute a buffer. The number of such buffer bands is governed by the input variable nbdbuf. It can happen that the highest (or two highest) band(s), if not separated by a gap from non-treated bands, can exhibit a very slow convergence rate. This buffer allows achieving convergence of important, non-buffer bands. In the present case, 6 bands have been converged with a residual better than tolwfr, while the two upper bands are less converged (still sufficiently for graphical representation of the band structure). In order to achieve the same convergence for all 8 bands, it is advised to use nband=10 (that is, 8 + 2). ## Using AbiPy to automate the most boring steps¶ The AbiPy package provides several tools to facilitate the preparation of band structure calculations and the analysis of the output results. First of all, one can use the abistruct script with the kpath command to determine a high-symmetry k-path from any file containing structural information (abinit input file, netcdf output files etc.). The high-symmetry k-path follows the conventions described in [Setyawan2010]. Let’s try with: abistruct.py kpath tbase3_5.in # Abinit Structure natom 2 ntypat 1 typat 1 1 znucl 14 xred 0.0000000000 0.0000000000 0.0000000000 0.2500000000 0.2500000000 0.2500000000 acell 1.0 1.0 1.0 rprim 0.0000000000 5.1080000000 5.1080000000 5.1080000000 0.0000000000 5.1080000000 5.1080000000 5.1080000000 0.0000000000 # K-path in reduced coordinates: # tolwfr 1e-20 iscf -2 getden ?? ndivsm 10 kptopt -11 kptbounds +0.00000 +0.00000 +0.00000 #$\Gamma$+0.50000 +0.00000 +0.50000 # X +0.50000 +0.25000 +0.75000 # W +0.37500 +0.37500 +0.75000 # K +0.00000 +0.00000 +0.00000 #$\Gamma\$
+0.50000  +0.50000  +0.50000 # L
+0.62500  +0.25000  +0.62500 # U
+0.50000  +0.25000  +0.75000 # W
+0.50000  +0.50000  +0.50000 # L
+0.37500  +0.37500  +0.75000 # K
+0.62500  +0.25000  +0.62500 # U
+0.50000  +0.00000  +0.50000 # X


To visualize the band structure stored in the GSR.nc file, use the abiopen script and the command line:

abiopen.py tbase3_5o_DS2_GSR.nc --expose -sns=talk


It is also possible to compare multiple GSR files with the abicomp script and the syntax

abicomp.py gsr tbase3_5o_DS1_GSR.nc tbase3_5o_DS2_GSR.nc -e -sns=talk


to produce the following figures:

For further details about the AbiPy API and the GSR file, please consult the GsrFile notebook .